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Abstract 

Nonlinear asymptotic integrators are applied to one-dimensional, nonlinear, autonomous, 
dissipative, ordinary differential equations. These integrators, including a one-step explicit, 
a one-step implicit, and a one- and a two-step midpoint algorithm, are designed to follow 
the asymptotic behavior of a system approaching a steady state. The methods require 
that the differential equation be written in an particular asymptotic form. This is always 
possible for a one-dimensional equation with a globally asymptotic steady state. In this 
case, conditions are obtained to guarantee that the implicit algorithms are well defined. 
Further conditions are deter min ed for the implicit methods to be contractive. These meth- 
ods are all first order accurate, while under certain conditions the midpoint algorithms may 
also become second order accurate. The stability of each method is investigated and an 
estimate of the local error is provided. 


1 Introduction 

Modeling non-equilibrium material behavior, such as viscoplasticity or chemical reactions, in- 
volves dissipation according to the second law of thermodynamics. The processes of primary 
interest are modeled by dissipative, ordinary differential equations with at least one asymptot- 
ically stable steady state. Most classical algorithms for the numerical solution of systems of 
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autonomous, ordinary differential equations are inefficient when used to predict the response of 
a material body as it approaches a steady state. The small step sizes that are usually required 
drastically increase the computational expense and time required to produce accurate solutions. 

Nonlinear asymptotic algorithms, which are intended to be more efficient in tracking asymp- 
totic behavior, can be defined by several types of asymptotic approximation. These algorithms 
were originally proposed for viscoplastic evolution problems by Walker and Freed (1991); 
such evolution is driven by thermodynamic dissipation. The primary goal of this paper is to 
investigate the consistency and stability of these algorithms and to determine when those that 
are implicitly defined have a unique solution. 

Dissipative evolution equations define a non-Hamiltonian dynamical system. In a Hamilto- 
nian system, the volume in state (phase) space is conserved over time. A dissipative system is 
one in which the volume in state space, as opposed to the material, contracts. The instantaneous 
rate of change of the local volume is the divergence of the velocity field. This could vary from 
point to point and may even be positive. Lichtenberg and Lieberman (1992) call a system 
dissipative if the local volume averaged over an orbit in state space contracts. In continuum 
mechanics, the condition that the trace of the Jacobian matrix of X = F(X) is negative implies 
that the density of the continua is increasing and, thus, that the volume is decreasing over time. 
Therefore solutions must behave asymptotically. 

A solution algorithm is most effective when it reflects the behavior of the original differential 
equation. The asymptotic algorithms are particularly appropriate for dissipative systems. A 
one-dimensional system, x = /(x), is dissipative if df /dx < 0 for all x. If such a system has a 
unique steady state, it is globally asymptotically stable; the exact solutions to the differential 
equation are asymptotic to the steady state. The asymptotic algorithms of Walker and Freed 
(1991) require that the differential equation be written in a special form. This is always possible 
for dissipative systems with a globally asymptotic steady state. Most of the results presented 
in this paper are for such equations. 

Two major issues are considered for the asymptotic algorithms applied to this family. First, 
the method must be well defined. In other words, there must exist a unique solution to the 
difference equation. If the algorithm is implicit, there are two numerical methods in common 
usage to solve the nonlinear difference equation for the value at the next time step, either the 
Newton-Raphson method or the method of successive approximations. The Newton- Raphson 
method, when applied to solve these equations, must approach the desired root of the difference 
equation. In general, restrictions on the step size are required to guarantee such convergence, 
see for example the Kantorovich theorem (Rall 1974). Furthermore, if the equation f(x) = 0 
has several roots, it is difficult to guarantee that the Newton-Raphson method will select the 
correct root. For example, for the equation x 3 — x = 0, perturbations of a starting value near the 
inflection point can cause the Newton-Raphson method to converge to any of the three roots. 
One solution is to apply the algorithms only in situations which guarantee that the implicit 
difference equation has a unique root. 

Second, any useful algorithm must converge to the exact solution for the given nonlinear 
differential equation. Convergence to a limit is controlled by the stability of the algorithm. 
The limi t, is the exact solution if the algorithm is first-order accurate. These two criteria are 
verified under some limits on the step size for the asymptotic algorithms applied to the famil y 
of dissipative differential equations. The generalized midpoint asymptotic algorithms discussed 
are shown to be second order accurate under certain conditions. 

The classic ideas of stability, A-stability, L-stability, B-stability, and G-stability, all were 
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defined in terms of systems of linear equations. Stability is largely determined by the location 
of the eigenvalues of the Jacobian of the linearized system in the complex plane. There is still 
no general agreement on a proper definition of stability for algorithms nonlinear in the step size 
when applied to nonlinear systems of ordinary differential equations. Some possible definitions 
axe considered here. 

The asymptotic algorithms are defined in Section 2. Dissipative, one-dimensional, ordinary 
differential equations are defined in Section 3, along with the special formulation required to write 
the difference equation for the asymptotic integrators. The conditions that these algorithms 
should satisfy are summarized in Section 4. The final four sections develop the properties of the 
asymptotic explicit, implicit, two-step midpoint, and the one-step midpoint algorithms. 


2 The Asymptotic Algorithms 

An asymptotic algorithm is intended to accurately model the behavior of a solution as it ap- 
proaches a steady state. In one dimension, then, the algorithms can be used to integrate 
autonomous differential equations x = f(x) such that the equation, f(x) = 0, has at least one 
root that is an asymptotically stable steady-state for the evolving system. 

The asymptotic integrators described in WALKER and FREED (1991) require that the one- 
dimensional, ordinary differential equation be written in the form 

x = c(x)[a(x) — x]. (1) 


Walker and Freed derived a variety of asymptotic approximations by first applying an inte- 
grating factor to Eqn. (1). 


rt+h 

x(t + h) = exp - J c(x(r])) drj 


x(t) 


exp -J c(x(f))df c{x(r]))a(x(r})) dr]. 


( 2 ) 


The first-order asymptotic integrators are obtained by approximating the integrals in Eqn. 
(2). To shorten the notation, put c„ = c(x n ), etc., and similarly for a(x), where x n is the 
approximation to the solution at the n time step, x(t n ). In particular, one obtains the first- 
order, asymptotic, forward integrator, by approximating c(x) by the constant and a(x) by 
a n . 

X n +l = X n + (l - e -Cnfc ) (an - X n ) , (3) 

and the first-order, asymptotic, backward integrator, by approximating c(x) by the constant c„ + i 
and a(x) by a n+ j. 

x»+i = x n + (l - e -c “+ lfc ) (a„ + i - x n ) . (4) 

A 0- midpoint integrator — first presented in Freed et al. (1992) — is equivalent to approxi- 
mating c(x) by the constant and a(x) by a n +* , resulting in the first-order , two-step , asymp- 
totic , 6 -midpoint integrator, 

Zn+1 = Xn + (l - e~ Cn+ ° h ') (an+e - x n ) , (5) 


3 



where 


x n +6 = X n + [1 - exp (-dcn+etyKcin+e - x n ), 
with Cn+e = c(x n+ e) and a n+e = a{x n +e). 

A one-step 0-midpoint integrator is developed below. While the two-step 0-midpoint inte- 
grator can be viewed as first applying the implicit and then the explicit asymptotic algorithm, 
the one-step method is obtained by first applying the explicit and then the implicit asymptotic 
algorithm. It is the first-order, one-step, asymptotic, 0-midpoint integrator, 

x n +i = x„ + (l - e -K 1 -*) c " + * c "+ 1 ]'‘) (a n - x n ) + (l - e ~ e ^ +lh ) (o^i - a n ). (6) 

Whenever 0 = 0, these reduce to the forward integrator, Eqn. (3); likewise, whenever 0 = 1, 
they reduce to the backward integrator, Eqn. (4). 

These asymptotic algorithms each produce the exact solution for a one-dimensional, linear, 
ordinary differential equation, x = c(a — x), where a and c are constants. The exact solution is 

x(t) = [x(0) — o]e _c< + a, 

while the algorithms yield 

x n+ i = [x n - a]e~ ch + a. 

On the other hand, the classical one-step forward and backward Euler methods do not yield the 
exact solution in this simplest of all situations. Therefore, the asymptotic algorithms can be 
expected to be an improvement over linear algorithms. 


3 Formulation of the System of Equations 

For one-dimensional equations, the definition of the algorithms requires that the equation be 
written in terms of an asymptote, a(x), and a coefficient, c(x). Conditions are required for the 
autonomous differential equation to be rewritten in the form of Eqn (1), for those cases when 
it is not already presented in this form. Irrespective of the form of a one-dimensional, ordinary 
differential equation, it can only have monotonic solutions. 

Lemma 3.1 All solutions to an autonomous, single-dimensional, differential equation, dx/dt = 
f(x) are monotonic. 

Proof: If any non-constant solution x(t) is not monotonic, then there is a t* such that x(t*) = x* 
for some x* such that x = /(x*) = 0. A root, x*, of f(x) = 0 produces the constant solution 
x(t) = x* when the initial condition is x(0) = x*. These two solutions passing through (f*,x*) 
violate the existence and uniqueness theorem for differential equations. 

This reflects the fact that a one-dimensional system can only have stable or unstable, isolated, 
steady states. Periodic orbits axe possible in two or higher dimensions, and chaotic behavior is 
possible in three or higher dimensions. 

A solution, x(t), approaches steady state if lim t _,oo dx/dt = 0. All such solutions to dx/dt = 
f(x) are asymptotic to a constant solution x(t) = x*, where /(x*) = 0. The sign of df /dx 
determines which constant solutions are asymptotes for other solutions, i. e. those which axe 
locally asymptotically stable. Only if df/dx < 0 at x* satisfying /(x*) = 0 is the root x* a 
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locally asymptotic steady state. In the formulation, dx/dt = f{x) = c(x)[a(x) - x], for c(x) > 0, 
the term a(x), the asymptote function, determines how many and which steady states exist, 
i. e. where a(x) - x = 0. For example, if a(x) = 0, then there is a unique steady state, and it is 
globally asymptotically stable. 

Lambert (1991) makes the following definitions. 

Definition 3.2 The system dx/dt = f(x) is dissipative in a region R x [a, b] if the norm 

(/(«) - f(y),x- y) < 0 


for all x,y G R and all t G [a, 6] • 

Definition 3.3 The system dx/dt = /(x) is contractive in the interval [a, 6] if two solutions 
x(t) and y(t ) with different initial conditions satisfy |x(i 2 ) — 3 /(^ 2 )! < l x (^i) — 3/(*i)l f or ^ 1^2 
satisfying a <t\ <t 2 <b. 

Any dissipative system is contractive because the difference |x(t) — y(t)\ is a non-increasing 
function of time when (/(x) - f(y), x-y) < 0. Consequently, any one-dimensional system with 
a globally asymptotic steady-state is contractive. 

Lemma 3.4 If the one- dimensional system dx j dt = /(x) has a globally asymptotic steady state, 
it is contractive. 

Proof: dx/dt = f(x) has a globally asymptotic steady state if df /dx < 0 for all x. Therefore 
x > y implies /(x) < f(y). In the one-dimensional case, the norm is 

[/(*) - /(y)K x - y) < °- 

Therefore the system is dissipative and also contractive. 

This result suggests that to reflect the behavior of the differential equation, the integrator 
should also be contractive. This idea will be discussed in general in the next section, and for 
the specific algorithms in later sections. 

The expression of a given system of ordinary differential equations in the form of Eqn. (1) 
is not likely to be unique. The two extremes of the range of equivalent expressions are obtained 
from putting all the nonlinearity in the coefficient c(x) by making the asymptote zero, or from 
putting all the nonlinearity in the asymptote by making the coefficient c(x) a constant. Numer- 
ical experimentation suggests that the more efficient technique is to put all the nonlinearity in 
the coefficient c(x), whenever possible. 

In a one-dimensional system, the existence of a globally asymptotic steady state is equivalent 
to being able to put all the nonlinearity in the coefficient, c(x). In higher dimensional systems, 
being able to put all the nonlinearity in the coefficients, c(x), implies that the system has a 
unique globally asymptotic steady state, but the converse is not true. If necessary, a coordinate 
change guarantees that x = 0 is the stable state. 

Lemma 3.5 Let dx/dt = f(x) be such that x = 0 is the unique root of f(x) and df /dx < 0 for 
all x, so that x = 0 is a globally asymptotic steady-state. Then the differential equation can be 
written in the form dx/dt = c(x)(0 - x), where a(x) = 0. Conversely, if dx/dt = -c(x)x, for 
c(x) > 0, then x = 0 is a globally asymptotically stable steady-state. 
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Proof: Define the coefficient, c{x), by 



if x / 0 

if x = 0 


Then c(x ) is continuous by L’Hospital’s rule. It is positive for all x since the fact that x = 0 is 
globally asymptotic implies that f(x) > 0 for x < 0 and f(x) < 0 for x > 0. 

Conversely, the Lyapunov function V (x) = x 2 shows that x = 0 is a globally asymptotic 
steady-state because dx/dt = 0 at x = 0 and dV/dt = 2 x(dx/dt) = -2c(x)x 2 < 0 when c(x) > 0. 


Example 3.6 That the formulation is not unique can he shown by the differential equation, 
dx/dt = —x 3 — x, which has a globally asymptotically stable equilibrium at x = 0. It is in the 
form of (1) with c(x) = 1 and a(x) = — x 3 . It can also can be written as dx/dt = — (x 2 + l)x 
so that c(x) = x 2 + 1 > 0, and a(x) = 0. Notice that the unique minimum of c(x ) occurs at the 
steady state in this example. 

Example 3.7 The equation x = 1 — x 3 has a globally asymptotic steady state at x = 1. The 
coordinate transformation y = x — 1 changes the equation into one having a unique steady state 
at y = 0, with dy/dt = —y 3 — y 2 — y. The coefficient, c(x), can be taken as either a constant, 
c = 1, in which case a(y) = — y 3 — y 2 , or with all the nonlinearity in the coefficient of y. In the 
latter case, c{y) = y 2 + y + 1 > 0, which reaches its absolute minimum at y = —1/2, a value 
unequal to the steady state, y = 0. 

Example 3.8 The requirement that the system have a globally asymptotic steady state in order 
to be written in the desired form cannot be dropped. The equation x = 6e -cx , for b and c positive 
constants, is dissipative but has no steady state. It cannot be rewritten in the form of Eqn. (1). 

The formulation of an autonomous, first-order, differential equation in terms of coefficients, 
c(x), and asymptotes, a(x), is not unique. Some judgement and analysis as to which formulation 
leads to the most efficient numerical modeling is required. 


4 Convergence of Nonlinear Algorithms 

To be useful in numerically integrating evolution equations, integration algorithms must satisfy 
some basic properties. If the algorithm is implicit, the algorithms must be shown to be well 
defined. The algorithm must also produce the exact solution in the limi t if the step sizes are 
taken very small. This convergence is usually a consequence of consistency and stability. An 
algorithm is stable if the approximation errors do not accumulate in such a way as to drive the 
approximate solution far from the exact solution after many steps. Consistency assures that an 
algorithm, which converges to some limit, actually converges to the exact solution. Finally, the 
algorithm should behave in the same way as the ordinary differential equation it is integrating. 
For example, if the differential equation has a globally asymptotic steady state, the algorithm 
should reproduce such behavior. It is in this property that the asymptotic integators distinguish 
themselves from other integrators for dissipative differential equations. 
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For an algorithm to be well defined, it must be determined that the function F{h, x n ) — £ n +i 
is an algorithm in the sense of ROSINGER (1980). When h is fixed, write Fh(x n ) = F(h,x n ). 
Let T be the largest value that the time step, h, can take. A compact subset, K, of the real 
numbers is a finite union of closed intervals. 

Definition 4.1 An algorithm is a family of mappings 

F h :X^X with he[0,T] 


on a topological space, X , such that 

1) F 0 = identity on X 

2) The function [0,T]xX-»X defined by (h,x) -*■ F h (x) is continuous in h and x. 

3) For each compact set, K C X, there is a M(K) > 0, such that for any pair u,v G K and 
h G (0, T], the local Lipschitz condition holds: 

|F fc (u) - F h {v)\ < [1 + M(K)h]\u - u|. 


Furthermore, the stability and consistency of the algorithms must be investigated to verify 
that the algorithm converges to the exact solution. An algorithm is convergent if the global 
discretization error, E = x(t + h) — x n +i tends to zero as h tends to zero. In other words, an 
algorithm converges if, as the step size tends to zero, the approximate solution approaches the 
exact solution. Consistency is defined for an arbitrary algorithm involving n degrees of freedom, 
X, by Marsden (1992). 

Definition 4.2 An algorithm Fh : 3f? n - ► 9?", for X = /(X), such that Fh(X n ) = X n+ i for 
h > 0 is called consistent or first order accurate if 


<mXn) 

dh 


= f(*n). 

h = 0 


Higher order accuracy is defined similarly by requiring that the higher order derivatives of 
Fh are equal to the corresponding derivative of x(t) as h tends to zero. 

This definition agrees with the more traditional statement of consistency in which the algo- 
rithm is written in the form 

x n +i = x n + h(f>{x ny K) 

and an algorithm is consistent if 

h — >-0 

The equivalence is seen by writing the Taylor series expansion for small h, 


F h (X n ) = X n + 


dF h (X n ) 

dh 


h + . . . . 

h=0 


The final requirement is that the algorithm be stable. Stability for linear algorithms can be 
investigated by examining the position of the eigenvalues in the complex plane. Classically, the 
stability of methods such as the various Runge-Kutta methods have been studied by applying 
the method to the one-dimensional test equation, x = Xx. A-stability, for example, is defined 
in terms of this linear test equation, see QuiNNEY (1987, p. 77). 
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Definition 4.3 A difference method is A-stable if when applied to the test problem, dx/dt = Ax, 
where Re{ A) < 0, all difference solutions, x n , tend to zero as n tends to infinity. 

The asymptotic algorithms produce the exact solution to this equation and are therefore more 
stable than A-stable. It might be more important to test the nonlinear asymptotic algorithms 
on nonlinear, ordinary differential equations, such as a dissipative system. 

Butcher proposed a definition of stability in 1975 to extend A-stability to Runge-Kutta 
methods applied to nonlinear, ordinary differential equations (Hairer and Wanner 1991, p. 
192). This is now called B-stability and can be applied to arbitrary methods, which may be 
nonlinear in the time step, h. 

Definition 4.4 A method is B-stable if when applied to contractive, ordinary differential equa- 
tions x = f(t,x) with (f(t,x) - f{t,y),x - y) < 0, then for all h> 0, |x n+ i - y„+i| < |x n - y n \. 

B-stability implies A-stability. For example, the linear implicit (backward) Euler method is 
B-stable. 

Lambert (1991) suggests that a useful definition for the stability of a nonlinear algorithm 
is that it be contractive, a generalization of B-stability to arbitrary algorithms. This condition 
reflects, in the algorithm, the fact that \y(t) — x(t)\ is a decreasing function of t if the differential 
equation is dissipative. 

Definition 4.5 If {x n } and {y n } are solutions generated from different initial conditions by an 
algorithm, then the algorithm is contractive if 

K+I -y n +i\ < \x n -y n \- 

Some authors, such as Hughes (1983), have defined stability to mean that |x n+ i| < |x„|. 
This condition is weaker than contractivity because it is implied by contractivity but does not 
imply contractivity. 

Lemma 4.6 Suppose an ordinary differential equation has a steady state x* = 0. If the algo- 
rithm is contractive, then |x n+ i| < |x n |. 

Proof: In |x n+i - y n+1 \ < |x n - y«|, let y n = y n+1 = 0. 

These ideas are generalized further in RosiNGER’s (1980) definition of stability for nonlinear 
algorithms. 

Definition 4.7 A difference scheme is stable if for each compact set, K C X, there is a L(K) > 
0, such that for any pair u,v € K, h E (0,T],n E N (where N is the set of positive integers), 
nh < T, the uniform local Lipschitz condition holds: 

\F£(u) - FZ(v)\ < L(K)\u - v\, 

where Ff is the n th iterate of Fh . 

Definition 4.8 The scheme is unconditionally stable if the algorithm is stable for h E (0, oo) 
and conditionally stable if it is stable only for h E (0,T], where T < oo. 

With these definitions, RosiNGER (1980) proved that for consistent algorithms, stability is 
equivalent to convergence. This is a generalization of the Lax-Richtmyer theorem for linear 
algorithms. 

Theorem 4.9 If an algorithm is first order accurate, then stability in the sense of Definition 
(4-7) is equivalent to convergence. 
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5 Explicit Algorithm 

The explicit asymptotic algorithm is clearly well defined. HENRICI (1962, p. 71) gives a theorem 
for the convergence of explicit methods in which a Lipschitz condition plays the role of stability. 

Theorem 5.1 Let the differential equation dx/dt = f(t,x ) have an approximate method defined 
by x n+ i = Zn + h<j>(t,x n ,h), where <f>(t,x,h ) is continuous in the domain D : t € [a,6],x € 
(- 00 , 00 ), and h € [0,/i o ], for h 0 > 0. Suppose there is a constant L such that \ <f>(t,x,h) - 
<f>(t, y, h)| < L\x—y\ for all (t, x, h) and ( t , y, h) in D. Then the consistency condition <p(t, x, 0) = 
f(t, x) is a necessary and sufficient condition for the convergence of the method defined by the 
increment function (f>. 

This theorem can be applied to the explicit, one-step, asymptotic algorithm for the au- 
tonomous system, x = /(x) = c(x)[a(x) — x]. Define </>(x,h) on the domain D : t € [a, b],x 6 
(- 00 , 00 ), and h € [0, h a ], for some h Q > 0 by putting 

4>{x,h) = {[a(x) - x][l - exp (-c(x)h)]}/h h € (0, h 0 ] 

4>{x,0) = f(x) h- 0. 

Lemma 5.2 Suppose dx/dt = f(x) = c(x)[o(x) - x], where c(x) and a(x) are continously 
differentiable for all x. Then <p(x,h) is continuous and continuously differentiable with respect 
to x on D. 

Proof: <p(x, h ) is clearly continuous in x for h > 0. By L’Hospital’s rule, 

lim <p(x, h ) = lim{[a(x) - x]c(x)e -c(x)fc } = /(x) = <p(x , 0). 

h—+0 h—> 0 

Therefore <f>(x, h ) is continuous in x, and h on D. 

Since a(x) and c(x) are assumed to be continuously differentiable for all x, for all fixed 
h € (0,hj, dfjdx exists and is continuous. Define df/dx = dffdx at each (x,0). Then by 
L’Hospital’s rule, for fixed x, at h = 0, 

lim d<t>/dx = Um{[a(x) — x](dc/dx)he~ c ^ h + [da/dx — 1][1 — e~ c ^ h )}/h 

= [a(x) - x](dc/dx) + [da/dx - l]c(x) = df /dx. 

Theorem 5.3 Suppose dx/dt = f{x) = c(x)[a(x) - x], where c(x) and a(x) are continously 
differentiable for all x. The explicit algorithm x n +i = x n + [a(x n ) — x n ][l — exp(— c(x n )h)] is 
convergent if \ <f>(x, h) — <f>(y, h)| < L|x — y| for all (x, h) and (y, h ) in D. 

Proof: The Henrici theorem gives the result immediately by the lemma, and because consistency 
holds, /(x) = <f>(x , 0). The algorithm is convergent without any restriction on the sign of c(x). 

In general, the method is not unconditionally convergent since h a must be chosen to make 
the Lipschitz condition in x hold for all (x, h) and (y, h) in D. This condition can be met if the 
derivative df /dx is bounded. 

The continuous function, 4>(x,h), satisfies the Lipschitz condition in x if df/dx is bounded 
for all x e (- 00 , 00 ), by the mean value theorem. In this case, the algorithm is convergent. The 
goal then is to choose h a so that 

d(f>/dx = {(a(x) - x]( dc/dx)he~* x)h + [da/dx - 1][1 - 
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is bounded for all a; € (— 00 , 00 ). 

This condition for convergence is not always satisfied, even if a global solution exists for 
dx/dt = f(x). For example, no such h a exists for either of the systems 

dx/dt = c(bx 2 — x) 

dx/dt = bx(a — x ), for a, b, c constants. 

The function, <f, associated with these differential equations does not satisfy the Lipschitz con- 
dition for all x € (— 00 , 00 ). 

5.1 Stability 

Dahlquist (1963) has shown that no explicit, linear, multistep method is A-stable. However, 
the explicit asymptotic algorithm is an exact solution for the test equation, and therefore, 

Lemma 5.4 The explicit, one-step, asymptotic algorithm is A-stable. 

Proof: For the equation dx/dt = Ax, the algorithm becomes 

Xn+l = x n e Xh , 

so that 

= W |e Ak | < |x„| |e*W‘|. 

Therefore 

|awi| < |Xo| |e* e < A) *| n+1 . 

But |e^ A)fc | < 1 for Re(X)h < 0, and so lim n _oc |x n | = 0. The method is A-stable. 

The stability properties of the explicit asymptotic algorithm axe therefore an improvement 
over any linear explicit method, such as the forward Euler. The explicit asymptotic algorithm is 
convergent when applied to nonlinear, ordinary differential equations satisfying a boundedness 
condition. 

Theorem 5.5 Suppose dx/dt = f(x) = c(x)[a(x)— x], where c(x) > 0 anda(x) are continuously 
differentiable for all x. The explicit algorithm x„+i = x n + [a(x n ) — x n ][l — exp(— c{x n )h)\ is con- 
vergent if df/dx is bounded for all x € (— 00 , 00 ). In this case, the algorithm is unconditionally 
stable. 

Proof: Since [1 — < c(x), and < 1, if c(x) > 0 and h > 0, 

d<f>/dx = {[a(x) — x](dc/dx)he~^ h + [(da/dx) — 1][1 — e~ c ^ h ]}/h 
< [a(x) — x]{dc/dx) + \{da/dx) — l]c(x) = df/dx. 

Therefore if df/dx is bounded for all x € (— 00 , 00 ), so is d<f>/dx. Then by the mean value theo- 
rem, the Lipschitz condition holds for <f>. By the Henrici theorem, the algorithm is convergent. 
It is unconditionally stable because there is no restriction on h. 
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Example 5.6 The equation, x = -x 3 - x, is a dissipative equation that is unbounded on the 
real numbers. If c(x) = 1 and a(x) = -x 3 , the explicit algorithm is both unstable and oscillates 
around the globally asymptotic stable-state for a step size ofh = — ln(0.5) = 0.693 and an initial 
condition greater than 1. On the other hand, if all the nonlinearity is collected in the coefficient 
c(x) = x 2 + 1 so that a(x) = 0, the explicit algorithm is stable in the sense that |x n+ i| < |x„| 
and does not oscillate about the steady state, for all h. 

In general, then, the algorithm is not stable for nonlinear, ordinary differential equations unless 
the step size is restricted. The stability of the explicit algorithm applied to a particular ordinary 
differential equation also depends on the formulation of that equation. 


5.2 Error Estimate 

The accuracy of the explicit asymptotic algorithm can be compared to the classical linear algo- 
rithm for small time steps. Estimates of the local truncation error can be obtained in the case 
that the ordinary differential equation dx/dt = f(x) is written in the form 

dx/dt — — c(x)x, 

with c(x) > 0. Its exact solution is from Eqn. (2), 

I rt+h 

- J c ( x {v)) dr) • 

The explicit, one-step, asymptotic algorithm, for this case, is x n +i = x n exp[-c(x n )/i]. Assume 
that x n = x(i). 

The truncation error is the difference of the exact and approximate solutions, 

E - x n+ i - x(t + h). 


It is more effective to first compute the ratio, 

[ [t+h 

x(t + h)/x n + 1 = exp - J c(x(r])) drj + c{x n )h 

By the mean value theorem, for some t < s < rj, 

c(x(rj)) - c(x n ) = ( dc(s)/drj ) (r) - t). 

Then 

I rt+h I ft+h 

jf c(x(T])) dr) - c{x n )h = jf [c{x(r])) - c{x n )] dr) 

I rt+h i 

J t ( dc{s)/dr )) (r) - t ) dr) 


< Mh 2 / 2, 


(7) 


where M is the maximum of \dc/dt\ on the closed interval [t,t + h]. Therefore, since e x > e~ 
the ratio is 

x(t + h)/x n+ 1 > exp(— Mh 2 /2). 
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then 


E < ®n+i [1 - exp(-Mh 2 /2)]. 

The local truncation error is proportional to h 2 for small h, and the cumulative error is propor- 
tional to h. This is the same error obtained for the classical forward Euler method. 

The explicit algorithm was called underdamped by WALKER and Freed (1991) in its ap- 
plication to a preditor-prey model. The algorithmic response depends on whether the functions 
c(x) and the exact solution x(t) are increasing or decreasing over the time interval of interest. 
For e xam ple, by Eqn. (7), the ratio, of the exact solution and the approximation is 


x(t + h)/x n+ 1 = exp 



- c(x n )]dr]\ . 


If x n > 0, x(r ) ) is decreasing and c(x) is increasing, then the difference c(x(r?)) - c(x n ) is non- 
positive for all -q. This implies that x(t + h)/x n+1 > 1. In this case, therefore, the explicit 
method approximation x n+ i is less than x(t + h). 


Example 5.7 The ordinary differential equation, x = -x 3 -x can be written asx = (x 2 +l)(-x) 
with c(x) = x 2 + 1. For the initial condition x n = 1 and h = 0.1, the exact solution is x(0.1) = 
0.832523, and the explicit asymptotic algorithm produces x n +i = 0.809675 . 


6 Implicit Algorithms 

The implicit solution, x n+ i, for any nonlinear implicit difference equation either may not exist or 
may not be unique. A numerical solution is usually sought at each step by either the Newton- 
Raphson method or by the method of successive approximations. Frequently, the Newton- 
Raphson method is the more efficient of the two methods. For implicit functions that have 
more than one solution, there is no guarantee to which root the Newton- Raphson method will 
converge, especially at larger time steps, so that it may not select the root closest to the exact 
solution to the differential equation. 

As an example in which the backwards Euler algorithm is not well defined, consider dx/dt = 
x 3 — x. This equation has constant solutions x = — l,x = 0, and x = 1, which are the only 
possible horizontal asymptotes. The graphs of the families of solutions can be sketched by 
looking at the signs of the first derivative, dx/dt = /(x) = x 3 - x, and of the second derivative, 
c ff/dt 2 = (df / dx)(dx / dt) = (3x 2 - l)(x 3 - x). The concavity of the solutions changes whenever 
x = ±l/-v/3. Of the possible steady states, the only one that is an asymptote is x = 0 because 
df fdx < 0 there. 

If |x*| < 1, so that it is in the region where solutions are asymptotic to x = 0, then the 
backwards Euler method is not well defined for any step size h > 0. In this case, x n+ i is a 
solution of the equation x n + fix 3 - (1 + h)x = 0. The discriminant of the cubic is negative if 
x 2 < 4(1 + h) 3 /27h. The right hand side takes its minimum of 1 at h = 1/2. Therefore there 
are three unequal real solutions to the backwards Euler difference equation for all h > 0 and 

|®n| ^ 1* 

For dx/dt = x 3 — x, in the form c[a(x) — x], so that c = 1 and a(x) = x 3 , the asymptotic 
algorithm gives 

Zn+1 = x n exp(— h) + X 3 +1 [1 - exp(-h)]. 
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The discriminant for this cubic is negative if x\ < 4/[27 exp(— h)(l — exp(— h))]. The right hand 
side takes its minim um value of 1 when h = ln(3/2). Therefore again, if |x n | < 1 there are 
three real and unequal solutions to the implicit asymptotic integrator difference equation for all 
h> 0. 

One might try to repair this situation by letting x„+i = F(h, x n ) be the root chosen by the 
Newton-Raphson method. In other words, the algorithm would be a combination of the differ- 
ence equation and the Newton-Raphson method. Unfortunately this function, x n+i = F(h,x n ), 
may not be continuous for fixed h and certainly will not always converge to an acceptable ap- 
proximation of the solution to the differential equation. The Newton-Raphson method does not 
always choose the root closest to the initial point. For example, let /(x) = x 3 — x = 0. Then 
if one starts at 0.5, the method converges to the root -1. Furthermore, the Newton-Raphson 
method behaves very poorly with perturbations of the initial guess. If one makes an initial guess 
0.4476, then the method converges to 1, if 0.4470 to 0, and if 0.4478 to —1. The Kantorovich 
theorem says only that Newton’s method converges to a root within the radius of convergence; 
it does not say the root is the only or desired one in that set. 

This problem can be avoided if the implicit algorithm has a unique solution. One condition 
which produces a unique solution is that the original differential equation have a globally asymp- 
totic, stable steady-state, x*. Assume that x n ^ x*, since if x n = x*, then the unique solution 
to the differential equation is known to be x(f) = x*, and no numerical method is required. This 
result is also valid for the classical backward Euler method. 

Lemma 6.1 Ifx = f(x) is dissipative and if f(x) = 0 has a unique solution, then the backward 
Euler method has a unique solution for all step sizes h. 

Proof: Define the function G(x) = x n + hf(x) - x. The backward Euler method has a unique 
solution if G(x) has a unique root. G(x) is monotonically decreasing because dG/dx < 0 if 
df jdx < 0. Assume that the unique steady state is x* = 0. The unique root of G(x) lies 
between x n and x* because G(x n ) = hf(x n ) and G(x*) = x n have opposite signs as a result of 
x* being a globally asymptotic steady state. Therefore G(x) has a unique root. 

In the case of a dissipative, ordinary differential equation with a unique steady-state, the im- 
plicit asymptotic algorithm has the same property as the differential equation. It has a globally 
asymptotic stable steady-state as h tends to infinity. Since lim/i_ 00 x n+ i = a(lim^_ 00 x m+1 ), and 
because x* is the unique value of x such that x* = a(x*), it is also true that lim/ l _ 00 x n+ i = x*. 


6.1 Existence and Uniqueness 

To verify that the asymptotic implicit algorithm is well defined for dissipative ordinary differ- 
ential equations, it must be shown that x n+ i exists and is unique. Define the function 

G{x) = x n e~* x)h + o(x)(l - e"^) - x (8) 

for fixed x n and h. The implicit difference equation has a unique solution for given x n and h iff 
G(x) has a unique root. The idea is to show that there is such a root between x n and x*. 

Le mma 6.2 Ifdx/dt = c(x)[a(x) — x] with c(x) > 0 has a globally asymptotic equilibrium state, 
x*, there is at least one root, r, of G(x) between x* and x n . 
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Proof: Either c(x*) = 0 or a{x*) = x* (or both). If c(x*) = 0, then 

G(x*) = x n — x*. 


If a(x*) = x*, then 

G(x*) = [x n - x*]e-« x "> h . 

Notice that the values of G(x*) agree if x* is a root of both c(x) and a(x) — x. Also, 

G(x n ) = [a(x n ) - x„][l - e -c(:,:n)fc ]. 

If x n > x*, then for the curve x(t) passing through x n to be asymptotic to x(t) = x*, dxjdt at 
x n must be negative and a(x n ) — x n < 0. Therefore G(x*) and G(x n ) have opposite signs, and 
by the continuity of G(x), a root must exist. A similar calculation produces the same result if 
x n < x*. 

Any root of G(x) must lie between x* and x„. 

Lemma 6.3 Suppose r is a root for G(x) and x n x* , then 

x* < r < x n if x* < x n 
x n < r < x* if x* > x n 

Proof: The root r ^ x* since G(x*) / 0. Because G(r) = 0, multiplying G(r) by c(r) and 
rearranging produces 

c(r)(a(r) — r) = c(r)[a(r) — x n ]e -< * r ^. 

Therefore dxjdt evaluated at r can also be written as 

dx/dt = c(r)[a(r) — x n ]e~^ h . 

If x* < r, then because x* is an asymptote and any solution must be decreasing, dxjdt evaluated 
at r is negative. The constraint c(x) > 0 implies that a(r) — x n < 0. Since G(x) can be rewritten 
as 

G(x) = [a(x) - x n ](l - e _c(x)fc ) - (x - x n ), 

G(r) = 0 implies that [a(r) — x n ](l _ e _c ( r ^) = r — x n , so that a(r) — x n and r — x n have the 
same sign. Therefore r < x n and x* < r < x n . The second statement follows in a similar 
manner. 

These two lemmas hold for x n in any region in which the solutions are asymptotic to a 
constant solution both from above and from below. They will be useful in proving results that 
do not hold for all x n and h. 

The root can now be shown to be unique by showing that G(x) is continuous and decreasing 
for all x near x*. This result is more restrictive than necessary for a specific x n because the 
monotonicity of G(x) is not required outside of the region between x n and x*. But it is required 
to account for arbitrary x n . The derivative of G(x), 

dGjdx = [a(z) — x n ](dc j dx)he~^ h + (dajdx)( 1 — e~^ h ) — 1 (9) 

is negative when h = 0, so that for G(x) to be monotonic, it must always be negative. The 
monotonicity of G(x) for all h can be easily verified if c(x) is a constant. 
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Theorem 6.4 Let dx/dt = f(x) = c[a(x)-x], where c > 0 is a constant, be such that f(x) = 0 
defines a unique globally asymptotic steady-state at the origin. Then the one-step, implicit, 
asymptotic algorithm has a unique solution for all step sizes h and all initial conditions. 


Proof: Because the origin is a globally asymptotic steady-state, df jdx < 0 for all x and conse- 
quently da/dx - 1 < 0. Since da/dx - 1 < 0 and since 0 < 1 - e -<A < 1, then 

dG/dx = (da/dx)( 1 - e -cfc ) - 1 < 0. 

Therefore the function G{x) is monotonic decreasing and has a unique root. 

If c(x) is not constant, then the solution may not be unique for all step sizes h. The behavior 
depends on the sign and magnitude of dc/dx. The following theorem covers the case when the 
coefficient c(x) is increasing for x > x* and decreasing otherwise. 

Theorem 6.5 Suppose that the ordinary differential equation dx / dt = c(x)(a(x) — x) withc(x) > 
0 , da/dx < 0 has a globally asymptotic, stable steady-state, x*, and that dc/dx and a(x)-x n have 
opposite signs when x is between x* and x n . Then the one- step, implicit, asymptotic algorithm 
has a unique solution for all time steps, h. 

Proof: This follows from the two previous Lemmas and the monotonicity of G(x). 

Example 6.6 The implicit algorithm applied to the differential equation, x = -x 3 - x written 
as x = — (x 2 + l)x, with c(x) = x 2 + 1 and a(x) = 0, has a unique solution due to this theorem. 

The behavior in general can be visualized from the fact that the term he vanes from 0 
when h = 0 to a maximum of e~ 1 /c(x) at h = l/c(x) and then asymptotically back to zero as h 
tends to infinity, for fixed x. Because (da/dx)( 1 - e~ c(x)h ) - 1 is negative for all h, in the case 
that [a(x) - x n ](dc/dx) is positive, there is a range of h near zero, and another range of large h, 
such that dG/dx < 0, so that the implicit algorithm has a unique solution for those h. 


6.2 Existence and Uniqueness by Contractions 

The use of the Newton-Raphson method to solve implicit nonlinear algorithms can lead to 
numerical difficulties. If the solution to the implicit algorithm can be shown to be the fixed 
point of a contractive function, then the method of successive approximations can often be used 
more efficiently in place of the Newton-Raphson method. The idea of a contractive mapping 
from a metric space to itself produces conditions under which an implicit algorithm can be 
shown to have a unique solution. 

De fini tion 6.7 A function f : M — *• M, where M is a complete metric space, is a contraction 

ff 

l/(x) - f(y) | < k\x - y\ for 0 < k < 1 

for all x,y € M. 
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A contraction has a unique fixed point, and the fixed point can be found by successive 
approximations starting at any initial point. 

An implicit algorithm is defined by the relation x n+ i = p(x n+1 ,x n ). For a given x n , the 
solution x n+ i is a fixed point of the function g(x). If the original differential equation has a 
globally asymptotic steady-state, then x n+ i must lie between x m and 0. Think of <?(x) defined on 
the complete metric space [0, x n ]. To prove that g{x) is a contraction on [0, x n ], it is sufficient to 
prove that the absolute value of dg/dx on [0, x n ] is less than one, since the mean value theorem 
says that |$(x) - g(y)\ = z ) |x - y\ for some x <z <y. 

For example, if all the nonlinearity is in the coefficient, c(x), then the implicit, one-step, 
asymptotic algorithm is x n +i = x n exp[— c(x n+ i)h]. In this case, g(x) = x n exp[— c(x)h]. It is 
then possible to restrict the time step, h, so that g(x) is a contraction. Let c* be the maximum 
of \dc/dx\ on the interval [0, x„]; c* exists because the interval is closed. 

\dg/dc\ = |x n (<2c/dx)hexp[— c(x)h]| < \x n (dc/ dx)h\ < |x n c*h| < 1 
Therefore the function g(x) is contractive on the interval if 

h < l/(x n c*). 

The more slowly that c(x) changes, the larger the step size possible. Under this restriction on 
h, the implicit algorithm has a unique solution, x„+i, the fixed point of g(x) and this point 
can be computed by successive approximations starting at, say, x„. This technique may be 
more efficient than the Newton-Raphson method, plus it guarantees that the solution obtained 
is unique. Notice further, that as the computed value of the solution gets closer to the steady 
state, 0, the size of the allowable step size increases. 


6.3 Stability 

Conditions have been found for the implicit difference equation to have a unique root, x n+1 , so 
that the algorithmic function F(h,x n ) = x n +i is well defined. Recall that Fh(x n ) = F(h,x n ) = 
x n+ i. To investigate the stability, one needs to know how the algorithm applied to dissipative 
differential equations behaves if the initial condition, x n , is perturbed to an initial condition, y n . 
The following discussion depends on the fact that x n can be written explicitly in terms of x n +i 
for fixed h as Fjf 1 (x n+ i) = x n , where, writing x for x n+ i, 

F^(x) = xe« x)h - afcXe*** - 1). 


Then 


dF^/dx = e* x)h + x{dc/dx)he* x)h 

~{da/dx)(e^ x)h - 1) - a(x)(dc/dx)he« x 
= e^ h - ( da/dx)(e^ h - 1) + [x - a(x)](dc/dx)he c ^ h . 


The following lemma will be useful in the investigation of stability. 


Lemma 6.8 Let dx/dt = /(x) = c(x)[a(x) — x], where c(x) > 0 and a(x) is monotonically 
decreasing, be such that /(x) = 0 has a unique globally asymptotic steady-state at the origin. 
Suppose either a) c(x) is a constant or b) dc/dx and a(x) — x n have opposite signs when x is 
between x* and x n . The unique root, x n+ x is an increasing function of x n for fixed h. 
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Proof: In case a) dcjdx = 0 implies that dF^/dx > 0. In case b), because da/dx < 0 and 
[x - a(x)](dc/dx) > 0 ,dF^/dx > 0. Therefore F^(x) is increasing and has an increasing 
inverse. This implies that x n+ \ = F^ (x n ) is an increasing function of x n . 

Lemma 6.9 Suppose a(x) is monotonically decreasing and either a) c(x ) is a constant or b) 
dcjdx and a(x) - x n have opposite signs when x is between x* and x n . Then the algorithm is 
contractive, |x n — y n \ > |x n +i — yn+i|> f or h. 

Proof: Suppose that x n > y n . Put x = x n +i and y = y n +i- Recall that x n = [x— a(x)]{exp[c(x)h] — 
1} + x. Then 

x n -yn = [x- a(x)]{exp[c(x)h] - 1} - [y - a(y)]{exp[c(y)h] - 1} + x - y. 

Verify that the function H(x) = [x — a(x)]{exp[c(x)/i] — 1} is increasing by showing that 

dH/dx = [x - a(x)](dc/dx)h exp[c(x)h] + [1 - (da/dx)]{exp[c(x)h] - 1} > 0. 

In case a), da/dx < 0 implies that dH/dx = [1 - (da/dx)] {exp[c(x)h] - 1} > 0. In case b), 
da/dx < 0 and [x - a(x)](dc/dx) > 0 imply that dH/dx > 0. Now x n > y n implies that x > y 
since x n +i is an increasing function of x n by the previous Lemma. H(x) increasing gives 

[x - a(x)]{exp[c(x)h] - 1} - [y - a(i/)]{exp[c(y)h] - 1} > 0. 

Therefore, x n - y n > x n+1 - y n+i > 0 and \x n - y n \ > \x n+x - y n+ 1 |, for all h. 

Any contractive algorithm satisfies Rosinger’s local Lipschitz condition. 

Theorem 6.10 Under the conditions that the implicit algorithm is contractive, it is also con- 
vergent. 

Proof: The result follows from RosiNGER’s theorem (4.9). The implicit algorithm is first-order 
accurate, or consistent, for any ordinary differential equation in the form x — c(x)[a.(x) x]. 

The derivative of x n+x with respect to h is 

dx n+x /dh = [—(dc/dx n+ i)(dx n +i/ dh)h - c(x n+ i)]{x n exp[-c(x n+ i)h] - a(x n+ i)} 

+ {1 - exp[-c(x n+1 )/i]}(da/dx n+ i)(dx n+ i/h). 

Therefore, the algorithm is first-order accurate because 

lim dx n+x /dh = c(x n )[a(x n ) - x n ] = x(x n ). 

h-+o 

Consequently, the algorithm is convergent whenever is it unconditionally stable. 


6.4 Error Estimate 

Estimates of the local truncation error for small step size h can be obtained in the case that the 
ordinary differential equation dx/dt = f(x) is written in the form 

dx/dt = — c(x) x, 
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with c(x) > 0. Its exact solution is 


/ r t+h 

c(x(rj))drj 


The implicit, one-step, asymptotic algorithm, for this case, is x n+ i = x n exp[— c(x n+ i)h]. Assume 
that x n = x(t). The ratio of the exact and approximate solutions is 

[ rt+h 

- j c(x(r])) drj + c(x n+ i)h . 

The local truncation error, E , is defined by x(t + h) + E = x n+1 . Since c(x) is continuous and 
differentiable, there is a bound L such that 

|c(x(t + h)) - c(x n+ i)| < L\x(t + h) - x n+ i| = L\E\. 

Then 

1 ft+h | ft+h 

\J c(x(tj)) dr] - c{x n+ i)h = jf [c(x(r})) - c(x(t + h)) + c(x(t + h)) - c(x B+ i)] dr) 

k t+h I ft+h 

(dc(s) /dr)) {t + h-rj)dr) + 1^ [c(x(t + h)) - c(x n+ i)] dr? 

< Mh 2 /2 + L\E\h, 

where M is the maximum of |dc/dt| on the closed interval [t,t 4- ft]. Therefore, the ratio is 

x(t -I- h)/x n+ 1 > exp(—Mh 2 /2 — L\E\h). 

Since E — 2 - 11 + 1(1 x(t -I- j x n ^.\^ 

E < x n+ i [1 — exp(— Mh 2 /2 — L\E\h)]. 

E is proportional to x n+1 (Mh 2 /2 + L\E\h). Therefore, the local truncation error is proportional 
to h 2 for small h, and the cumulative error is proportional to h. This is the same error obtained 
for the classical backward Euler method and for the first-order, forward, asymptotic algorithm. 

The implicit, asymptotic algorithm was said to produce overdamping by Walker and 
Freed (1991). 

Lemma 6.11 If x = — c(x)x and if c(x) > 0 is increasing, then whenever x n > 0, x n +i > 
x(t + h). 

Proof: Assume the contrary, that x n+ i < x(t + h). Since x n > 0, the exact solution is decreasing. 
In the integral, 

[ ft+h 

-J [c( x (v)) - c(x n+ i)]drj , 

the term 

c(x(??)) - c(x n+ i) = [c(x(r?)) - c(x(t + h ))] + [c(x(t -I- h )) - c(x n+ i)] 
is positive since c(x) is increasing. Therefore 

x(t + h ) /x n+1 < 1 and x(t + h) < x n +i, 
which is a contradiction. 

Example 6.12 The ordinary differential equation, x = — x 3 — x can be written as x = (x 2 + 
1)(— x) with c(x) = x 2 + 1. For the initial condition x n = 1 and h = 0.1, the exact solution is 
x(0.1) = 0.832523, and the implicit asymptotic algorithm produces x n+ i = 0.842796 . 
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7 Two-Step Midpoint Algorithm 

The classical, linear, midpoint method approximates the integral of x = f(x) from x n to x n +i 
by a rectangle having base h and height taken at the midpoint of the interval [i, t 4- h]. HUGHES 
(1983, p. 145), generalized this to a nonlinear algorithm, in terms of x n and x n+1 , by taking the 
height at x n + a = (1 — a)x n + OLX n+ \. In other words, x n+a is a linear interpolation between x n 
and x n+ i. The nonlinear midpoint algorithm of HUGHES (1983, p. 144) and of HOGGE (1978) 
was applied to equations of the form, M(x,t)(dx/dt) + K(x,t)x = F(t), by substituting x n+Q 
to get M n+a v n+a + K n+a x n+a = F n+a , where v = dx/dt. Because x n+Q is a linear function of 
x n and x n +i, the expressions in terms of n + a drop out, and one solves for x n 4 -i. 

The asymptotic midpoint algorithm of Eqn. (5) can be viewed as a combination of the results 
of one single-step algorithm on part of the interval followed by another single-step algorithm on 
the remainder of the interval. If one uses the result from the implicit integrator, Eqn. (4), taken 
over the interval [t,t+6h] as the initial condition for the explicit integrator, Eqn. (3), taken 
over the remaining interval [t -\-9h, one obtains an algorithm in which both Xn+i and the 

0-midpoint, x n+e , must be determined. This is a two-step algorithm. The midpoint algorithm 
can be expressed as 

x n+ i = x n exp (~c(x n+e )h) + [1 - exp(-c(x n+ #)h)]a(x n+9 ), 

in agreement with Eqn. (5), where, for c^+e = c(x n+ g) and dn+e = a { x n+e), 

X n +e = x n exp (~dc n+e h) + [1 - exp(-0c n+e /i)]a n+e . 

The value x n+e exists under the same restrictions applying to the implicit asymptotic algo- 
rithm; it is obtained from the implicit method applied over the interval [t, t+6h]. In those cases, 
x n+1 exists since it is defined explicitly from x n+e . For example, if c(x) is a constant, then 

Zn+i = x n exp(-ch) + [1 - exp(-ch)]a(x n+ #), 

and the midpoint algorithm is well defined for all 0 < 9 < 1. 

The midpoint algorithm, when 0 = 1/2, is more accurate than any of the other two-step 
0-midpoint algorithms since, for small h , it approximates the second derivative of the solution 
as well as the first derivative. 

Lemma 7.1 The two-step midpoint algorithm is first-order accurate for all 0 < 0 < 1. It is 
second-order accurate iff 6 = 1/2. 

Proof: The derivative of x n+ i with respect to h is 

dx n+ i/dh = [(dcn+g / dx)xOh - c^+g] exp (-c„ + eh)(x n - o^+o) + [1 - exp(-<^ + *h)](da n+ */dx)x0. 


Letting h tend to zero, 


dx n +i 

dh 


h = 0 


CnX n + CjiOji — X n . 


Therefore the algorithm is first-order accurate for all 0 < 0 < 1. 
derivatives with respect to h shows that 


Likewise, talcing second 


d?x n +i 


dh 2 


h=0 


20(dc n /dx)(a n - x n )x n + 26c n (da n /dx)x n - c„x n 
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But 


X — [fi(f j dx)x — (dCnj X,f)X n "i" ^n^n) 

where x = f(x). Therefore 

(Px n+ i I 


dh 2 


= x n iff 0 = 1/2. 


h = 0 


Only when 0 = 1/2 is the midpoint algorithm second-order accurate. 


7.1 Stability 

This midpoint asymptotic algorithm reduces to the trapezoidal algorithm defined by Freed et 
al. (1992) if the function a(x) is identically zero. This is the form of the ordinary differential 
equation that HUGHES (1983) took as the test case. Hughes tests the stability of his algorithm 
by application to the model equation, 

dx/dt + c(x, t)x = 0 for c(x,t ) > 0. 

If this system is linear, his trapezoidal and midpoint algorithms agree. Stability is defined by 
Hughes as requiring |x n+ i| < |x n |. His midpoint algorithm is unconditionally stable if a > 1/2, 
and if a < 1/2, it is stable whenever h satisfies 


Cn+ih < 2/(1 - 2a). 


As Hughes points out, this shows that while the trapezoidal algorithm is unconditionally stable 
when applied to linear problems, it is not when applied to nonlinear problems. The asymptotic 
midpoint algorithm, on the contrary, is unconditionally stable for the same values of a when 
applied to both linear equations and to the nonlinear family with all nonlinearity in the coefficient 
of x. 

If the asymptotic midpoint algorithm is applied to Hughes’ model equation dxjdt + c(x)x = 0 
for c(x) > 0, which has a unique steady-state x = 0, then 


x n 4-i — x n exp[ c(x n .|_0)/i]. 


No matter what c(x n +e) is, |x n+ i| < |x n | for all h and 0. Therefore, in this weaker definition 
of stability, the asymptotic midpoint algorithm, in the case where the asymptote is zero, is 
unconditionally stable for all 0 < 0 < 1, at least under the Hughes definition of stability. As 
such, it is a more stable algorithm than the Hughes, nonlinear, midpoint algorithm. 

7.2 Local Truncation Error 

For small h , the linear midpoint algorithm has a local truncation error of order h 3 , while the for- 
ward Euler method has a truncation error proportional to h 2 for small h. A similar improvement 
is provided by the asymptotic midpoint algorithm. 

Lemma 7.2 The local truncation error for the midpoint algorithm applied to x = — c(x)x is 
proportional to h 3 iff 0 = 1/2. Otherwise, it is proportional to h 2 . 
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Proof: Assume that x(t) = x n . The ratio of the exact to the approximate solution over the 
interval [t , t + h] is 

[ rt+h 

- j c(x(rj)) dr] + c(x n+0 )h . 

Take the Taylor series of c(x(r])) about t 0 = t + Oh, 

c(r 7 ) = c(t 0 ) + (77 — t 0 )(dc(t 0 ) / dt) + [(77 — te) 2 /2\(d 2 c(£)/dt 2 ), 
for some ( such that t 0 < £ < 77 . The local truncation error requires the computation of 

rt+h. 

T n = J c(rj)dr] — c(x n+0 )h 

= c(t 0 )h - c(x n+ e)h + ^(77 - t 0 )(dc{t 0 )/dt)dr] + ^[(v - t 0 ) 2 /2}{d 2 c(t)/dt 2 )d V 
= c{te)h - c(x n+0 )h + [(1 - Of/ 2 - O 2 /2]h 2 (dc(t 0 )/dt) + 0(h 3 ). 

Since c(x) is continuous and differentiable, there is an L such that 

c(#n+ 0 ) c (^(^)) = L\x n +0 “ 


Therefore, 

T n = LE e h + [(1 - ef/2 - O 2 /2]h 2 c'(t 0 ) + 0(h 3 ), 

where E 0 is the order h 2 , local error from the asymptotic, implicit algorithm. T n is of order h 3 
if and only if the coefficient of h 2 is zero. 

(1 - Of /2 - 0 2 / 2 = 0 iff 0 = 1/2. 

Since the local truncation error is E < x n+x [l - exp(-|T n |)], E is proportional to h 3 iff 

0 = 1 / 2 . 

This agrees with the previous results that the truncation error for both the asymptotic 
explicit ( 0 = 0) and implicit (0 = 1) methods is of order h 2 . 


8 One-Step Midpoint Algorithm 

A one-step, generalized, midpoint method is defined in which first the explicit asymptotic algo- 
rithm and then the implicit method is applied. The order of application of the explicit and the 
implicit algorithms is the opposite of that employed to define the two-step, generalized, midpoint 
algorithm. The explicit asymptotic integrator, Eqn. (3), taken over the interval [t,t+ (1 - 0)h] 
gives the initial condition for the implicit asymptotic integrator, Eqn. (4), taken taken over the 
remaining interval [t+ (1 - 0)h,t+h], The form of the algorithm is 

Sn+l = x n + (1 - exp{ — [(1 - 0)c n + OCn+ilh}) (dn - X n ) + (1 - eXp[-0C„ +1 h]) (<Xn+l - On)- (10) 

When 0 = 0, the explicit asymptotic algorithm is recovered, and when 0 = 1, the implicit 
asymptotic algorithm is recovered. 
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If the time coefficient, c, is a constant, and the asymptote is monotonically decreasing, the 
one-step, asymptotic, midpoint algorithm is uniquely defined. If c is constant, the one-step 
midpoint algorithm can be written as 

x n+ i = (x n - a n ) exp(-ch) + o„ exp(-0ch) 

+ [1 - exp(— 0ch)]on + i. 

If in addition a is a constant so that the original equation is linear, 

z»+i = [*» “ a]e -di + a. 

The one-step midpoint algorithm gives the exact solution for linear equations, just as the other 
asymptotic integrators do. 

Theorem 8.1 If the dissipative equation x = c[a(x) — x] has c > 0 constant and dajdx < 0, 
then the one-step midpoint algorithm is well defined. 

Proof: Define the auxiliary function, 

G(x) = (x n — 0 ^) exp(— ch) + a* exp(— 8ch) 

+ [1 — exp(— 0ch)]a(x) — x. 

The derivative, dG/dx = [1 — exp (—0ch)](da/dx) — 1 < 0 for all h and 6 so that G is monotonic 
in x for any fixed h. If G(x) crosses the x axis, this root is unique. 

G(x n ) = [l cxp( c/i)](o n x n ). 

G(x*) = (x„ - c^) exp (-ch) + [a„ - a(x*)] exp (~6ch), 

since a(x*) = x*. Because x* is a globally asymptotic steady-state and because a(x) is decreasing, 
the sign of the terms, (x n — a n ) and [a n — a(x*)] are the same and opposite to (a n — x n ). Therefore 
G(x n ) and G(x*) have opposite signs, and the unique x n+ i lies between x n and x*. 

If the function a(x) is identically zero, this midpoint algorithm also reduces to the trapezoidal 
algorithm of Freed et al. (1992). In this case, 

x n -j-i = x n exp[ (1 O^Cnh Ocn+xH^. (11) 

Theorem 8.2 The one-step midpoint algorithm applied to the system x = — c(x)x, with c(x) > 
0, is well-defined for all h> 0 and 8. 

Proof: Put 

G(x) = x n exp[— (1 — 9)cnh] exp[— 9c(x)h] — x. 


Then for x* = 0, 


G{x n ) = -x n (l - exp[— c(x n )h]) 

G(x*) = x n exp[— (1 — ^JcnhJexpf— 0c(O)h] 



have opposite signs. Therefore if G(x) is monotonic between x n and 0, it must have a unique 
root and the algorithm is well defined. This calculation depends on the facts that since the 
differential equation has a globally asymptotic steady-state, 

— c(x) — ( dc/dx)x < 0; 


and 

9hexp[— 9c(x)h] < l/c(x). 

If x n > x > 0, then 

dG/dx = ® n exp[— (1 — 6)c(x n )h][—(dc/dx)6h]ex'p[—9c(x)h] — 1 

< [c(x)0h] exp[— (1 - 0)c(x n )h ] exp[-0c(®)h] - 1 

< exp[— (1 - 9)c(x n )h ] — 1 < 0. 

Therefore there is a unique value x n+i satisfying the asymptotic, one-step, midpoint algorithm 
for all h and 6. 

The accuracy of this midpoint algorithm is also improved when 6 = 1/2 since, for small h , 
it appro xim ates the second derivative of the solution as well as the first derivative. 

Lemma 8.3 The one-step midpoint algorithm is first-order accurate for all 0 < 0 < 1. It is 
second-order accurate iff 6 = 1/2. 

Proof: The derivative of x n+i with respect to h is, by Eqn. (10), 

^±i = [(1 - 0)c n + 0Cn + i + 0(dc/dx)(dx n+ i/dh)h] (1 - exp{-[(l - 0 )cn + (a n - x n ) 

dh 

+ [0c„+i -1- 0(dc/dx)(dx n+1 /dh)h] exp[-^c n+ ih](a„ + i - a„) 

+ (1 - expt-tfCu+i/i]) (da/dx)(dx n+ if dh) 


Letting h tend to zero, 


dx n+ i 

dh 


— Cn{&n ®n) — %n- 

h=0 


Therefore the algorithm is first order accurate for all 0 < 9 < 1. 
derivatives with respect to h shows that for arbitrary 0, 


Likewise, taking second 


(PXn+l 

dh 2 


— 20(dcn/dx)(a n x n )x n *b 2 9c^ n fda^ l jdx)x n c n x n 

h = o 


But 

where ® 


x — ( df/dx)x = (dc n /dx)(a n - x n )x n + Cnidan/dx^n - c»i n , 
f(x). Therefore 


d?x. 


n - hi 


dh 2 


= X* 


fc=0 


iff 9 = 1/2. 


When 0 = 1/2, the one-step midpoint algorithm is second-order accurate. 
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8.1 Stability 

If the asymptotic, one-step, midpoint algorithm is applied to dx/dt + c{x)x = 0 for c(x) > 0, 
which has unique steady-state x = 0, then x n+ i is given by Eqn. (11). If one puts Cn + e = 
(1 9)c(x n ) “I - $c(x n -n.), then 

X n +l = Zn exp[— Cn+fl/l] 

In any case, |x n+1 | < |x n | for all h and 9. Therefore, like the two-step method, the asymptotic, 
one-step, midpoint algorithm, in the case that the asymptote is zero, is unconditionally stable 
in this weak sense for all 0 < <j> < 1. 


8.2 Local Truncation Error 

The local truncation error, for small h, improves to order h 3 for the one-step, asymptotic, 
midpoint algorithm iff 0 = 1 / 2. 


Lemma 8.4 The local truncation error for the one-step, asymptotic, midpoint algorithm applied 
to x = —c(x)x is proportional to h 3 iff 9 = 1/2. Otherwise, it is proportional to h 2 . 


Proof: Assume that x(t) = x n . The ratio of the exact to the approximate solution over the 
interval [t, t + h] is 

I ft+k 

- ^ c(x(r})) drj + (1 - 9)c(x n )h + 0c(x n+ i)h . 

The integral in the exponent can be expanded using Taylor series. 

rt+h 


c (x( 7 ?)) dr) + (l- 9)c(x n )h + 0c(x n +i)h 

= / [c(*(77)) - Cn] drj + / [c(x(t?)) - 0c» +1 ] dr) 

Jt Jt+(l-6)h 

dc(t) [*+( 


J +( } [(v -t) + 0(2)1 dr) + / [(v -t-h) + 0(2)}dr] 

Jt dt Jt+ll-0)h 


dt 


^(1 - 9) 2 h 2 / 2 - dc{t + h) 9 2 h 2 /2 + 0(h 3 ) 


dc{t ) , 

dt 

dc(t) 

dt 


[(i - efh 2 ! 2 - e 2 h 2 / 2] + 


dc(t) dc(t + h) 


dt 


dt 


9 2 h 2 / 2 + 0(h 3 ) 


[(1 - 0) 2 h 2 j2 - 9 2 h 2 / 2] + 0(h% 


since dc(t)/dt is continuous and differentiable implies that dc(t)/dt — dc(t + h)/dt is proportional 
to h. Therefore, the integral is proportional to h 3 iff 6 = 1/2. 

Since the local truncation error is E < x n+1 [l — exp(— \T n \)], E is proportional to h 3 iff 

9 = 1 / 2 . 

This also agrees with the previous results that the truncation error for both the asymptotic 
explicit ( 6 = 0) and implicit (0 = 1) methods is of order h 2 . 
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9 Summary 

Three nonlinear asymptotic algorithms proposed by WALKER and FREED (1991) have been 
investigated for one-dimensional, dissipative, ordinary differential equations of the form, x = 
c(x)[a(x) — x] for c(x) > 0 with a globally asymptotically stable-steady state. A new one-step, 
midpoint, asymptotic algorithm has been proposed. 

A method related to those studied here was developed by Dennis (1959) to integrate ho- 
mogeneous, first-order, ordinary differential equations with exponential type solutions. In the 
equation x = /(x), the substitution u = ln(x) transforms the ordinary differential equation to 
u = f/x. The algorithm x r+p = x r exp[/ gdt], where g = f/x, is obtained in which a numerical 
solution is taken for the integral. If the exact solution is exponential, then u varies slowly, and 
the algorithm gives good results. 

To apply these asymptotic algorithms, the ordinary differential equation must be written in 
the form, x = c(x)[a(x)— x]. Such a formulation is not unique, and the choice affects the behavior 
of the algorithms. If the equation is dissipative and has a globally asymptotic steady-state, then 
the equation can always be written in the form, x = — c(x)x with c(x) > 0. 

The algorithms are particularly studied for two extreme cases of dissipative differential equa- 
tions with a globally asymptotic steady-state: that c(x) is constant or that a(x) — 0. The 
differential equation with c(x) a constant is a generalization of the semi-linear equations. The 
implicit algorithm is well-defined for any step size h if c(x) is constant. In the second case, a 
bound on h was found for the implicit integrator to have a unique solution and, in fact, to be 
contractive. In the first case, the implicit algorithm is contractive for all h > 0. They may be 
efficiently solved with the method of successive approximations. 

Each of the asymptotic algorithms produces the exact solution for any linear differential 
equation and are, therefore, more successful than the classical linear methods. In the case that 
a(x) = 0, both midpoint integrators, and the trapezoidal integrator proposed in Freed et al. 
(1992), reduce to the identical algorithm. 

All of the algorithms are first-order accurate. Both midpoint methods are second-order 
accurate iff 6 = 1/2. They use a combination of the explicit and implicit one-step integrators. 
However, they differ from a predictor-corrector method in that the component integrators are 
not both applied over the full interval h. The two-step midpoint algorithm has local truncation 
error of the order, h 3 , rather than the truncation error of order h 2 of the other asymptotic 
algorithms, but only if 6 = 1 /2. The one-step explicit and implicit asymptotic algorithms have 
the same order of local truncation error, h 2 , as the corresponding linear methods. 

Evaluation of the stability of nonlinear algorithms applied to nonlinear, ordinary differential 
equations is a difficult problem. Each of these asymptotic algorithms is A-stable. This in itself 
is an improvement over many linear algorithms. Conditions have been found under which the 
asymptotic algorithms are contractive, which is often taken as the definition of stability. The 
methods are not uniformly stable under the strongest definitions of stability. The range of 
values of the time step, h, that give stable behavior depends on the form in which the ordinary 
differential equation is written. In the case that a(x) = 0, they are uniformly stable in the sense 
that |x n -|_i| ^ 

The algorithms allow larger step sizes within the stable range near the steady state for 
differential equations that exhibit asymptotic behavior. As desired, the asymptotic algorithms 
behave comparably to a differential equation with a globally asymptotic steady-state. The 
algorithms will also be shown, in a subsequent paper, to be more than competitive with linear 
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multi-step methods when applied to stiff systems of equations. 
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